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Abstract: We present a fabrication process of low-cost superlattices and simulations related with the heat 
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PACS: 68.65.Cd, 66.70.Df, 81.16.-c, 65.80 -g, 31.12.xv 



Introduction 

Understanding and controlling the thermal properties of 
nanostructures and nanostructured materials are of 
great interest in a broad scope of contexts and applica- 
tions. Indeed, nanostructures and nanomaterials are get- 
ting more and more commonly used in various 
industrial sectors like cosmetics, aerospace, communica- 
tion and computer electronics. In addition to the asso- 
ciated technological problems, there are plenty of 
unresolved scientific issues that need to be properly 
addressed. As a matter of fact, the behaviour and relia- 
bility of these devices strongly depend on the way the 
system evacuates heat, as excessive temperatures or 
temperature gradients result in the failure of the system. 
This issue is crucial for thermoelectric energy-harvesting 
devices. Energy transport in micro and nanostructures 
generally differs significantly from the one in macro- 
structures, because the energy carriers are subjected to 
ballistic heat transfer instead of the classical Fourier's 
law, and quantum effects have to be taken into account. 
In particular, the correlation between grain boundaries, 
interfaces and surfaces and the thermal transport prop- 
erties is a key point to design materials with preferred 
thermal properties and systems with a controlled 
behaviour. 



In this article, the prediction tools used for studying 
heat transfer in low-cost superlattices for thermoelectric 
conversion are presented. The technology used in the 
fabrication of these superlattices is based on the method 
developed by Marty et al. [1,2] to manufacture deep sili- 
con trenches with submicron feature sizes (Figure 1). 
The height and periodicity of the wavelike shape of the 
surfaces can be monitored. When the trenches are filled 
in with another material, they give rise to superlattices 
with rough interfaces. This was the motivation for 
studying both the thermal conductivity and the Kapitza 
resistance [3] of superlattices with rough interfaces. We 
focus mostly at the influence of interfacial width of the 
superlattices made of two semiconductor-like materials, 
with simple Lennard-Jones potential for the description 
of interatomic forces. Simulations of the Kapitza resis- 
tance for binary system of silicon with metal are also 
presented. These interfaces are difficult to be modelled, 
first of all because of the phonon-electron coupling that 
occurs at these interfaces and secondly because of the 
plethora of potentials which can be used. The choice of 
potential is based in a comparison of their performance 
to predict in a correct manner, the harmonic and anhar- 
monic properties of the material. Results on the Kapitza 
resistance of a silver/silicon interfaces are also 
presented. 
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Figure 1 SEM pictures obtained by the group ESYCOM and ESIEE at Marne-la-Vallee, France, showing two submicron trenches in a 
silicon wafer 



Fabrication process of superlattices 

To reduce the processing time and the manufacturing 
costs, vertical build superlattices are proposed as 
opposed to conventional planar superlattices. In Figure 2, 
a schematic representation of the two types of superlat- 
tices is given comparing their geometries. With this pro- 
cess, silicon/metal superlattices can be fabricated. 
Although final device will have material layers in the 
tens of nanometre range, 5- and 15-um width superlat- 
tices are fabricated using typical UV lithography. These 
thick layer superlattices are necessary to develop an 
accurate model of thermal resistance at the metal/semi- 
conductor interfaces. 



Vertical superlattices were obtained by patterning and 
then etching the silicon by deep reactive ion etching 
(DRIE). The trenches were filled using electrodeposition 
on a thin metallic seed layer. In Figure 3, a scanning 
electron microscope (SEM) image of a processed silicon 
wafer with micro-superlattices is given. There are voids 
at the bottom of the trenches which are explained by 
the absence of the seed layer at the bottom, and the fact 
that they prevent any copper growth. These voids were 
successfully eliminated by increasing the amount of seed 
layer sputtered in subsequent trials. The excess copper 
on top, resulting from the trenches being shorted to 
facilitate electroplating, was polished away using 



Conventional Superlattices 



Vertical Superlattices 
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Planar deposits of SC 
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by MBE 
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by DRIE etching. 



Figure 2 Structural comparison between conventional superlattices and vertical superlattices. 
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Figure 3 SEM image of copper-filled 5-um-wide trenches 



chemical-mechanical polishing. This is done to electri- 
cally isolate the trenches from one another so as to 
allow thermo-electrical conversion. 

We aim to fabricate optimized vertical nano-superlat- 
tices (with layers ranging <100 nm each) with high ther- 
moelectric efficiency. High thermoelectric efficiency 
occurs for high electrical conductivity and low thermal 
conductivity. The electronic conductivity will be con- 
trolled though the Si doping and the use of metal to fill 
in the trenches. The film thickness needs to be 
decreased, to decrease the individual layer thermal con- 
ductivity and increase the influence of the interfacial 
thermal resistance. 

To obtain such dimension on a large area at low cost, 
we are developing a process based on the transfer by 
DRIE of 30-nm line patterns made of di-block copoly- 
mers [4]. For this purpose, it is required to characterize 
them to the best possible degree of accuracy. Measure- 
ments at this scale will possibly be plagued by quantum 
effects [5,6]. That is the reason why we fabricated first 
micro-scale superlattices, to make thermoelectric mea- 
surements free from quantum effects and then applied 
the method to characterize the final nano-superlattice 
thermoelectric devices. 

Simulations: thermal conductivity of superlattices 

When the layer thickness of the superlattices is compar- 
able to the phonon mean free path (PMFP), the heat 



transport remains no longer diffusive, but ballistic 
within the layers. Furthermore, decreasing the dimen- 
sions of a structure increases the effects of strong inho- 
mogeneity of the interfaces. Interfaces, atomically flat or 
rough, impact the selection rules, the phonon density of 
states and consequently the hierarchy or relative 
strengths of their interactions with phonons and elec- 
trons. Thus, it is important to study and predict the 
heat transfer and especially the influence of the height 
of superlattice's interfaces on the cross and in-plane 
thermal conductivities. This is a formidable task, from a 
theoretical point of view, as one needs to account for 
the ballistic motion of the phonons and their scattering 
at interfaces. Molecular dynamics is a relatively simple 
tool which accounts for these phenomena, and it has 
been applied successfully to predict heat-transfer prop- 
erties of superlattices. 

Two routes can be adopted to compute the thermal 
conductivity, namely, the non-equilibrium (NEMD) [7] 
and the equilibrium molecular dynamics (EMD) [8]. In 
this article, we have considered both methods to charac- 
terize the thermal anisotropy of the superlattices. In the 
widely used direct method (NEMD), the structure is 
coupled to a heat source and a heat sink, and the result- 
ing heat flux is measured to obtain the thermal conduc- 
tivity of the material [9,10]. Simulations are held for 
several systems of increasing size and finally thermal 
conductivity is extrapolated for a system of infinite size 
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[11,12]. The NEMD method is often the method of 
choice for studies of nanomaterials, while for bulk ther- 
mal conductivity, particularly that of high conductivity 
materials, the equilibrium method is typically preferred 
because of less severe size effects. Comparisons between 
the two methods have been done previously, concluding 
that the two methods can give consistent results [13,14]. 
Green-Kubo method for nanostructures is proven to 
have greater uncertainties than those of NEMD, but a 
correct description of thermal conductivity with EMD is 
achieved by establishing statistics from several results, 
starting from different initial conditions. 

The superlattice system under study is made of super- 
position of Lennard-Jones crystals and fee structures, 
oriented along the [001] direction. The molecular 
dynamics code LAMMPS [15-17] is used in all the 
NEMD and EMD simulations. The mass ratio of the 
two materials of the superlattice is taken as equal to 2, 
and this ratio reproduces approximately the same acous- 
tic impedance difference as that between Si and Ge. Per- 
iodic boundary conditions are used in all the three 
directions. Superlattices with period of 40a 0 are dis- 
cussed, where a 0 is the lattice constant. The shape of 
the roughness is chosen as a right isosceles triangle. The 
roughness height was varied from one atomic layer (1 
ML = l/2« 0 ) to 24«n- For each roughness, heat transfer 
simulations with NEMD were performed for several sys- 
tem sizes in the heat flux direction to extrapolate the 
thermal conductivity for a system of infinite size [11]. 
For EMD simulations, the size of the system is smaller 
than with NEMD simulations and only one size is con- 
sidered 20a 0 x 10a o x 40a 0 > where the last dimension is 
perpendicular to interfaces. 

In Figure 4, we gathered the results for the in-plane 
and cross-plane thermal conductivities obtained by the 
two methods. The thermal conductivity is measured 
here in Lennard-Jones units (LJU), which correspond in 
real units typically to W/mK. At the low temperatures 
considered (T = 0.15 LJU), the period of the superlattice 
is comparable to that of the PMFP. The qualitative 
interpretation of the results shows that the thermal con- 
tact resistance of the interface has a strong influence on 
the superlattice thermal conductivity. The results pre- 
viously obtained by NEMD method [12], and, in particu- 
lar, the existence of a minimum for the in-plane thermal 
conductivity are now confirmed using the EMD method. 
The evolution of the TC as a function of the interfacial 
roughness is found to be non-monotonous. When the 
roughness of the interfaces is smaller than the superlat- 
tice's period, the in-plane thermal conductivity first 
decreases with increasing roughness. It reaches a mini- 
mum value which is lower by 35-40% compared to the 
thermal conductivity of the superlattice with smooth 
interfaces. For larger roughness, the thermal 



conductivity increases. The initial decrease of the in- 
plane thermal conductivity is quite intuitive if one con- 
siders the behaviour of phonons at the interfaces, which 
may be described by two different models. In the acous- 
tic mismatch model [18,19], the energy carriers are 
modelled as waves propagating in continuous media, 
and phonons at the interfaces are either transmitted or 
specularly reflected. For atomically smooth interfaces, it 
is assumed that phonons experience mainly specular 
scattering. The roughness enhances diffuse scattering at 
the interface in all space direction. 

In the diffuse-mismatch model, on the other hand, 
phonons are diffusively scattered at interfaces, and their 
energy is redistributed in all the directions [20]. In prac- 
tice, the acoustic model describes the physics of interfa- 
cial heat transfer at low temperatures, for phonons 
having large wavelengths, while the diffuse model is 
relevant for small wavelengths phonons. At the consid- 
ered temperature in the current study, we are most 
probably in an intermediate situation where the physics 
is not captured by one single model. Nevertheless, both 
models predict that a moderate amount of interfacial 
roughness will tend to decrease the in-plane TC, 
because rough interfaces will increase specular reflection 
and diffusive scattering of phonons travelling in the in 
plane direction. However, if the roughness is large 
enough, then locally, the phonons encounter smooth- 
like interfaces, and the partial group of phonons that are 
diffusely scattered in all space direction decreases. This 
might explain the further increase of the thermal con- 
ductivity when the roughness is large enough. 

The behaviour of the cross-plane thermal conductivity 
is different: it increases monotonously with the interfa- 
cial roughness. For smooth interfaces, the cross-plane 
thermal conductivity is 50% lower than the in-plane 
thermal conductivity. This anisotropy has to be taken 
into account for thermal behaviour of systems made of 
sub-micronic solid layers. Invoking again the acoustic 
mismatch model, we conclude that the transmission 
coefficient of the solid/solid interface is smaller than the 
reflection coefficient, which is not surprising if we con- 
sider the acoustic impedance ratio of the two materials. 
Roughness increases the transmission coefficient as it 
increases the diffused scattering at the interface [12]. 

The same qualitative trend regarding the influence of 
the roughness on the thermal conductivity of superlat- 
tices has been reported previously for materials with dif- 
fusive behaviour, without thermal contact resistance 
[21]. In this case, the variation of the in-plane and 
cross-plane conductivities with the interfacial roughness 
is due to the heat flux line deviation that minimizes the 
heat flux path in the material that has the lower thermal 
conductivity. This tends to increase the cross-plane 
thermal conductivity. On the other hand, the increase of 
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Figure 4 Cross-plane and in-plane thermal conductivity functions of the height of interfaces calculated by EMD and NEMD methods. 



the roughness leads to the heat flux constrictions that 
decrease the in-plane thermal conductivity. The qualita- 
tive interpretation of the results shows that the thermal 
contact resistance of the interface has a strong influence 
on the superlattice thermal conductivity. 

Simulations: Kapitza resistance 
Superlattices with rough interfaces 

The discussion above shows that obviously the phononic 
nature of the energy carriers has to be taken into 
account to understand heat transfer in superlattices, and 
that the evolution of the superlattice TC may be qualita- 
tively understood in terms of interfacial or Kapitza resis- 
tance. At a more quantitative level, the Kapitza 
resistance is defined by 



Rk 



AT 

T 



and thus quantifies the temperature jump AT across 
an interface subject to a constant flowing heat flux /. In 
general, the Kapitza resistance may be computed using 
NEMD simulations by measuring the temperature jump 
across the considered interface. For superlattices, how- 
ever, the direct method can be used only to measure 
easily the Kapitza resistance only for smooth surfaces, 
because of the difficulty involved in measuring locally 
the temperature jump for non-planar interfaces. To 
compute the Kapitza resistance for superlattices with 



rough interfaces, we have used EMD simulations, and 
the relation between R K and the auto-correlation of the 
total flux q(t) flowing across an interface: 



1 



SkriT 2 



/ {mm) 



dt 



where S is the interface area. The latter formula 
expresses the fact that the resistance is controlled by the 
transmission of all the phonons travelling across the 
interface. 

In the situation of interest to us here, the transmission 
of phonons is expected to be strongly anisotropic, and 
thus the resistance developed by an interface should 
depend on the main direction of the heat flux. To mea- 
sure this anisotropy, we have generalised the previous 
equation and introduced the concept of directional 
resistance, by considering the heat flux q%{i) ln the 
direction 6 in (0,n/2) with the normal of the interface. 

The resistance in the direction 8 may be then quanti- 
fied by the generalised Kapitza resistance: 



Re Sk B T 2 



+0O 

J (<fo(t)<fe(0))dt 



This angular Kapitza resistance quantifies the trans- 
mission of the heat flux in the direction making an 
angle 9 with the normal of the interface. 
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Figure 5 displays the generalised Kapitza resistances 
measured with MD for superlattices with rough inter- 
faces having variable roughnesses. Again, the results are 
displayed in LJU, which correspond to a resistance of 
10 x 10" m K/W in SI units. The period of the super- 
lattice considered is larger than the PMFP, which here 
is estimated to be around 20a 0 . We have focused on 
two peculiar orientations 0 = 0 and 9 = tt/2 which cor- 
respond, respectively, to the cross-plane and in-plane 
directions of the superlattices. It is striking that, for a 
given interfacial roughness, the computed resistance 
depends on the orientation 6. We have found that for 
almost all the systems analysed, the Kapitza resistance is 
larger in the cross-plane direction than in the direction 
parallel to the interfaces. This is consistent with the 
observation that the thermal conductivity is the largest 
in the in-plane direction (Figure 4). Again, this rein- 
forces the message that the heat transfer properties of 
superlattices are explained by the phononic nature of 
the energy carriers, and that theses energy carriers feel 
less friction in the in-plane direction than that in the 
direction normal to the interfaces. Measuring the direc- 
tional Kapitza resistance is a first step towards a quanti- 
tative measurement of the transmission factor of 
phonons depending on their direction of propagation 
across an interface. 
Silver/silicon interfaces 

The Cu and Ag films on Si-oriented substrates are the 
principal combinations in large-scale integrate circuits. 
Furthermore, with the fabrication process of vertical- 



built superlattices described in previous section, we are 
interested in the heat transfer phenomena related to the 
metal/semiconductor interfaces. The prediction of heat 
transfer in these systems becomes challenging when the 
thickness of the layers reaches the same order of magni- 
tude as the PMFP. For heat transfer studies, MD is well 
suited for dielectrics since only phonons carry heat. For 
metals, coupling between phonons and electrons can be 
modelled with the two-temperature model [22]. For the 
above systems, it has been proven that the Kapitza resis- 
tance is mainly due to phonon energy transmission 
through the interfaces [23,24]. The interfacial thermal 
resistance, known as the Kapitza resistance [25,26] is 
important to be studied as it might become of the same 
order of magnitude than the film thermal resistance. In 
this section, interatomic potentials for Ag and Si are dis- 
cussed. Using NEMD simulations, for an average tem- 
perature of 300 K, the Kapitza resistance of Si/Ag 
systems is determined. 

Modified embedded-atom method (MEAM) is the only 
appropriate potential that can be used for metal/semi- 
conductor systems. The first nearest-neighbour MEAM 
(INN MEAM) potential by Baskes et al. [27] and the sec- 
ond nearest-neighbour MEAM (2NN MEAM) by Lee 
[28] are examined in the current study. The general 
MEAM potential is a good candidate for simulating the 
dynamics of a binary system with a single type of poten- 
tial. For example, it can be applied for both fee and bec 
structures. Furthermore, this potential includes direc- 
tional bonding, and thus can be applied for Si systems. 
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In dielectric materials heat transfer depends mainly on 
phonons' propagation and their interactions. To make 
the best choice among a great number of potentials for 
calculating thermal conductivity, the dispersion curves 
and the lattice expansion coefficient were studied. Elec- 
tron transport predominates at the heat transfer in 
metals. MD cannot simulate electron movement, 
although some models are suggested in the literature to 
include the interactions between electron and phonons 
but without yet a satisfying results for investigating heat 
transfer. As it is not possible to test the quality of elec- 
tronic interactions, only the lattice properties are com- 
mented to determine the correct potential for 
simulating Ag. The dispersion curves in the [cj, 0, 0], [%, 
cj, 0] and [cj, cj, i;] directions are determined and com- 
pared with the experimental dispersion curves of Ag 
[29] for the INN MEAM and 2NN MEAM (Figure 6). 

To compare the anharmonic properties of Ag, the 
equilibrium lattice parameter is simulated for different 
temperatures using the INN MEAM, and 2NN MEAM 
potentials. This is modelled with an fee slab consisting 



of 108 atoms of silver with periodic boundary conditions 
in all the directions. Initially, the temperature of the 
crystal was 0 K. For each temperature the simulations 
are performed with a 20 ps constant-pressure simulation 
(NPT) during which the volume of the box occupied by 
the atoms for each temperature is stored. The mean 
value of the volumes of the equilibrated energy is used 
to calculate the linear expansion coefficient. For each 
constant temperature, the volume of the simulation box 
is divided by the volume at 0 K. This ratio is directly 
proportional to the expansion coefficient. The expansion 
coefficients of Ag, obtained for the two potentials are 
compared to the experimental values [30] in Figure 7. 
The uncertainties on the linear expansion coefficient 
variation are less than 5% compared with the experi- 
mental values. 

The 2NN MEAM potential allows recovering the 
expansion coefficient for Ag quite accurately while the 
INN MEAM potential significantly underestimates it. 
For Ag, the two potentials provide a good description 
for the more basic properties, such as cohesive energy, 
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lattice parameters and bulk modulus [31]. Even if the 
INN MEAM potential gives results closer to the experi- 
mental values for dispersion curves, the values obtained 
for the linear thermal expansion are not reasonable. 
Therefore, the INN MEAM potential cannot be consid- 
ered appropriate for simulating heat transfer for silver. 
Regarding the investigation of heat-transfer temperature, 
the 2NN MEAM gives the best results for harmonic and 
anharmonic properties for silver and for silicon using 
the previous results of the literature [32]. Kapitza resis- 
tance is predicted for the 2NN MEAM Si/Ag potential. 
The interface thermal resistance, also known as Kapitza 
resistance, R K , creates a barrier to heat flux and leads to 
a discontinuous temperature, AT", drop across the 
interfaces. 

The interactions between silicon and silver are 
described thanks to the 2NN MEAM potential in which 
the set of parameters has been determined to produce a 
realistic atomic configuration of interfaces. The model 
structure consists of two slabs in contact: one of Si with 
a diamond structure, and one of Ag. The periodic 
boundary conditions are used in all the directions and 
the Si crystal is composed of 7200 atoms, while the Ag 
crystal is composed of 2560 atoms. In the first stage of 
MD simulation, the system is equilibrated at a constant 
temperature of 300 K for 20 ps using an integration 
time step of 5 fs. The heat sources are placed in the 
extremes of the structure, and one layer of Si and Ag is 
frozen to block the movement of Si atoms in the z- 



direction. The temperature gradient is formed in the z- 
direction, imposing hot and cold temperatures above 
and below the fixed atoms in z-direction. Using an inte- 
gration time step of 5 fs, the simulation is run for 
5.0 ns, with an average system temperature of 300 K. In 
Figure 8, the temperature profile for the Si/Ag system is 
shown. 

The Kapitza resistance obtained with NEMD is 4.9 x 
10" m K/W. The temperature profile for Si is almost 
flat due its high thermal conductivity. With MD simula- 
tions, it is not possible to simulate heat transfer due to 
the electrons, and thus the steep slope of Ag is due to 
its low lattice thermal conductivity. The value i? KX is in 
the range 1.4-125 x 10" 9 m 2 K/W which also includes 
the Kapitza conductance for dielectric/metal systems 
[33,34]. 

Conclusions - Discussion 

A new fabrication method for superlattices is used, 
reducing the time and fabrication costs. With the fabri- 
cation of vertical superlattices, several questions a rose 
for the influence of the roughness' height of the super- 
lattices and the quality of interface on the thermal trans- 
port. When the length of the superlattice's period is 
comparable to the phonon-free mean path, the heat 
transfer becomes ballistic. 

The cross-plane and in-plane thermal conductivities of 
a dielectric/dielectric (representing Si-Ge systems) 
superlattice are predicted using EMD and NEMD 
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simulations. Both methods give the same tendencies for 
the anisotropic heat transfer at superlattices with rough 
interfaces. The in-plane thermal conductivity exhibits a 
minimum for a certain interfacial width, while the 
cross-plane thermal conductivity increases modestly in 
increasing the width of the interfaces. The Kapitza resis- 
tance of these interfaces is also studied, with a proposed 
methodology in this article, introducing the concept of 
directional thermal resistance. Values presented here are 
coherent with the difference between the in-plane and 
cross-plane thermal conductivities. 

Molecular dynamics simulations are also used to study 
the metal/semiconductor interfaces. Among all the 
interatomic potentials that are available, the MEAM 
potential is a good alternative to work with since it can 
be used for different materials. At 300 K, the 2NN 
MEAM potential gives the best results for the funda- 
mental properties associated with the heat transfer of 
silicon and silver. Previous results [23,24,32] suggest 
that interfacial thermal conductance depends predomi- 
nantly on the phonon coupling between silicon and 
metal lattices so that Si/Ag can be simulated without 
considering the contribution of electron heat transfer. 
The value of magnitude of the Kapitza resistance for a 
Si/Ag system is within the range of Kapitza resistance 
proposed in the literature. 

This study proves that making rough instead of 
smooth interfaces in superlattices is a useful way to 
decrease the thermal conductivity and finally to design 



materials with desired thermal properties. Furthermore, 
when more interfaces are added (rough or smooth), i.e. 
when the superlattice's period decreases, the interfacial 
thermal resistance becomes comparable to the superlat- 
tice's layers thermal conductivity. With these two para- 
meters, namely, the introduction of rough interfaces and 
the decrease of the superlattice's period, we can create 
systems with controlled values of the thermal 
conductivity. 
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